******************
*** 1_PART ONE ***
******************
cd "~"
use 0_prepared.dta, clear

	// APPENDIX
	// Sample composition 
	keep if !missing(lrscale,gener5,sex,partisan,intr)
	drop leftright
	label drop leftright
	gen leftright = .
	replace leftright = 1 if lrscale >= 0 & lrscale <= 2
	replace leftright = 2 if lrscale >= 3 & lrscale <= 4
	replace leftright = 3 if lrscale == 5
	replace leftright = 4 if lrscale >= 6 & lrscale <= 7
	replace leftright = 5 if lrscale >= 8 & lrscale <= 10
	label define leftright 1 "Left" 2 "Center-left" 3 "Center" 4 "Center-right" 5 "Right"
	label values leftright leftright

	bysort year: tab leftright [aw=anweight]
	bysort year: tab gener5 [aw=anweight]
	bysort year: tab sex [aw=anweight]
	bysort year: tab partisan [aw=anweight]
	bysort year: tab intr [aw=anweight]

	** Figure A1a
	graph hbar [aw=anweight], over(leftright) over(year) stack asyvars percentage ///
		ylab(, nogrid angle(h)) ytick(0(10)100) ytitle("") scheme(s1mono) plotregion(margin(l=2 r=2))  ///
		legend(pos(6) rows(1) stack region(lcolor(none)) symplacement(center) placement(center) j(center)) ///
		bar(1, fcolor(edkblue) lcolor(white) lwidth(thin)) ///
		bar(2, fcolor(edkblue*.8) lcolor(white) lwidth(thin)) ///
		bar(3, fcolor(edkblue*.6) lcolor(white) lwidth(thin)) ///
		bar(4, fcolor(edkblue*.4) lcolor(white) lwidth(thin)) ///
		bar(5, fcolor(edkblue*.2) lcolor(white) lwidth(thin)) name(dist_leftright, replace)

	** Figure A1b
	graph hbar [aw=anweight], over(gener5) over(year) stack asyvars percentage ///
		ylab(, nogrid angle(h)) ytick(0(10)100) ytitle("") scheme(s1mono) plotregion(margin(l=2 r=2))  ///
		legend(pos(6) rows(1) stack region(lcolor(none)) symplacement(center) placement(center) j(center)) ///
		bar(1, fcolor(edkblue) lcolor(white) lwidth(thin)) ///
		bar(2, fcolor(edkblue*.8) lcolor(white) lwidth(thin)) ///
		bar(3, fcolor(edkblue*.6) lcolor(white) lwidth(thin)) ///
		bar(4, fcolor(edkblue*.4) lcolor(white) lwidth(thin)) ///
		bar(5, fcolor(edkblue*.2) lcolor(white) lwidth(thin)) name(dist_gener, replace)

	** Figure A1c
	graph hbar [aw=anweight], over(sex) over(year) stack asyvars percentage ///
		ylab(, nogrid angle(h)) ytick(0(10)100) ytitle("") scheme(s1mono) plotregion(margin(l=2 r=2))  ///
		legend(pos(6) rows(1) stack region(lcolor(none)) symplacement(center) placement(center) j(center)) ///
		bar(1, fcolor(edkblue) lcolor(white) lwidth(thin)) ///
		bar(2, fcolor(edkblue*.8) lcolor(white) lwidth(thin)) name(dist_sex, replace)

	** Figure A1d
	graph hbar [aw=anweight], over(partisan) over(year) stack asyvars percentage ///
		ylab(, nogrid angle(h)) ytick(0(10)100) ytitle("") scheme(s1mono) plotregion(margin(l=2 r=2))  ///
		legend(pos(6) rows(1) stack region(lcolor(none)) symplacement(center) placement(center) j(center)) ///
		bar(1, fcolor(edkblue) lcolor(white) lwidth(thin)) ///
		bar(2, fcolor(edkblue*.8) lcolor(white) lwidth(thin)) name(dist_partisan, replace)

	** Figure A1e
	graph hbar [aw=anweight], over(intr) over(year) stack asyvars percentage ///
		ylab(, nogrid angle(h)) ytick(0(10)100) ytitle("") scheme(s1mono) plotregion(margin(l=2 r=2))  ///
		legend(pos(6) rows(1) stack region(lcolor(none)) symplacement(center) placement(center) j(center)) ///
		bar(1, fcolor(edkblue) lcolor(white) lwidth(thin)) ///
		bar(2, fcolor(edkblue*.8) lcolor(white) lwidth(thin)) ///
		bar(3, fcolor(edkblue*.6) lcolor(white) lwidth(thin)) ///
		bar(4, fcolor(edkblue*.4) lcolor(white) lwidth(thin)) name(dist_intr, replace)

	** Table A1
	use 1_p1_cntry.dta, clear
	collapse env_im env_hom env_un im_hom im_un hom_un cc_env cc_im cc_hom cc_un egal_redis, by(land year)

	tabstat env_im env_hom env_un im_hom im_un hom_un cc_env cc_im cc_hom cc_un egal_redis, stats(mean N) by(year)

// Constraint over time across countries
use 1_p1_cntry.dta, clear

foreach var of varlist env_im env_hom env_un im_hom im_un hom_un cc_env cc_im cc_hom cc_un egal_redis {
	gen `var'_val = .
	replace `var'_val = abs(`var')
}

// Constraint
egen rm = rowmiss(env_im_val env_hom_val env_un_val im_hom_val im_un_val hom_un_val)
egen pcA = rowtotal(env_im_val env_hom_val env_un_val im_hom_val im_un_val hom_un_val) if rm == 0
gen BSA = pcA*(2/(4*(4-1)))
lab var BSA "Constraint within belief system {it:A}"
drop rm

egen rm = rowmiss(env_im_val env_hom_val env_un_val im_hom_val im_un_val hom_un_val cc_env_val cc_im_val cc_hom_val cc_un_val)
egen pcB = rowtotal(env_im_val env_hom_val env_un_val im_hom_val im_un_val hom_un_val cc_env_val cc_im_val cc_hom_val cc_un_val) if rm == 0
gen BSB = pcB*(2/(5*(5-1)))
lab var BSB "Constraint within belief system {it:B}"
drop rm

egen rm = rowmiss(egal_redis_val)
egen pcec = rowtotal(egal_redis_val) if rm == 0
gen BSec = pcec*(2/(2*(2-1)))
lab var BSec "Constraint within economic belief system"
drop rm

collapse BS*, by(land year)

save 1_part_one.dta, replace

	// APPENDIX
	// Country-specific levels of constraint
	** Figure A6
		clonevar land_nomiss = land if !inlist(land,1,23,26,31,39)
		twoway (line BSA year, lc(cranberry)) (line BSec year, lc(edkblue*.5) lp(shortdash)) (line BSB year, lc(edkblue) lw(medthick)), ///
			by(land_nomiss, style(compact) note(" ")) ylab(, labsize(medlarge)) xlab(2002(20)2022, angle(90) nogrid labsize(medlarge)) ///
			xmtick(2002 2016 2022) xtitle(" ") ///
			legend(pos(6) rows(1) label(1 "Belief System {it:A}") labe(2 "Economic Belief System") label(3 "Belief System {it:B}"))

levelsof year, local(years)

foreach var of varlist BS* {
	gen `var'_mean = .
	gen `var'_lci  = .
	gen `var'_uci  = .

	foreach i of local years {
		capture {
	ci means `var' if year == `i'
		replace `var'_mean = r(mean) if year == `i'
		replace `var'_lci = r(lb) if year == `i'
		replace `var'_uci  = r(ub) if year == `i'
	}
}
}

bysort land: tab year if BSA != .
bysort land: tab year if BSB != .
bysort land: tab year if BSec != .
	
** Figure 1
twoway (bar BSB_mean year, color(ltblue*.7) lc(ltblue) barwidth(1.5)) (rcap BSB_lci BSB_uci year, lc(edkblue) lwidth(vthin) msize(vlarge)) ///
	(bar BSA_mean year, color(red*.5) lc(red)) (rcap BSA_lci BSA_uci year, lc(cranberry) lwidth(vthin) msize(medsmall)), ///
	xlab(2002(20)2022, nogrid labsize(10-pt)) xmlab(2004(2)2020, labsize(8-pt)) ylab(, labsize(10-pt)) xtitle(" ") ytitle("") ///
	legend(keygap(4) order(- "Belief System {it:A}:" 3 "Mean" 4 "95% CI" - "Belief System {it:B}:" 1 "Mean" 2 "95% CI") ///
	pos(6) rows(2) symplacement(center) symxsize(5) forcesize size(10-pt)) 

tabstat BSA BSB BSec if year == 2002, stat(mean sd)
tabstat BSA BSB BSec if year == 2012, stat(mean sd)
tabstat BSA BSB BSec if year == 2016, stat(mean sd)
tabstat BSA BSB BSec if year == 2022, stat(mean sd)

gen year0212=0 if year == 2012
	replace year0212 = 1 if year == 2002	
gen year1222=0 if year == 2022
	replace year1222 = 1 if year == 2012
gen year0222=0 if year == 2022
	replace year0222 = 1 if year == 2002
gen year1622=0 if year == 2022
	replace year1622 = 1 if year == 2016

** Table 3
ttest BSA, by(year0212)
ttest BSA, by(year1222)
ttest BSA, by(year0222)
ttest BSA, by(year1622)
ttest BSB, by(year1622)
ttest BSec, by(year0212)
ttest BSec, by(year1222)
ttest BSec, by(year0222)
ttest BSec, by(year1622)

ttest BSB = BSA if year == 2016
ttest BSB = BSA if year == 2020
ttest BSB = BSA if year == 2022

bysort year: ttest BSec = BSA
bysort year: ttest BSB = BSA
bysort year: ttest BSB = BSec

// Preparation for Figures 2a-2d
use 1_part_one.dta, clear

levelsof land, local(countries)

foreach var of varlist BSA BSB BSec {
	gen `var'_b = .
	gen `var'_se  = .
	gen `var'_ll  = .
	gen `var'_ul  = .
	gen `var'_n  = .
	gen `var'_beta = .

	foreach l of local countries {
			quietly {
				capture {
		reg `var' year if land == `l', beta
		matrix beta = e(beta)
			replace `var'_b = _b[year] if land == `l'
			replace `var'_se = _se[year] if land == `l'
			replace `var'_ll = _r_lb[year] if land == `l'
			replace `var'_ul = _r_ub[year] if land == `l'
			replace `var'_n = e(N) if land == `l'
			replace `var'_beta = beta[1,1] if land == `l'
			}
		}
	}
}

gen BSA_b_0212 = .
gen BSA_se_0212 = .
gen BSA_n_0212 = .
gen BSA_beta_0212 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSA year if inrange(year,2002,2012) & land == `l', beta
	matrix beta = e(beta)
		replace BSA_b_0212 = _b[year] if land == `l'
		replace BSA_se_0212 = _se[year] if land == `l'
		replace BSA_n_0212 = e(N) if land == `l'
		replace BSA_beta_0212 = beta[1,1] if land == `l'
		}
	}
}

gen BSA_b_1222 = .
gen BSA_se_1222 = .
gen BSA_n_1222 = .
gen BSA_beta_1222 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSA year if inrange(year,2012,2022) & land == `l', beta
	matrix beta = e(beta)
		replace BSA_b_1222 = _b[year] if land == `l'
		replace BSA_se_1222 = _se[year] if land == `l'
		replace BSA_n_1222 = e(N) if land == `l'
		replace BSA_beta_1222 = beta[1,1] if land == `l'
		}
	}
}

gen BSA_b_1622 = .
gen BSA_se_1622 = .
gen BSA_n_1622 = .
gen BSA_beta_1622 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSA year if inlist(year,2016,2020,2022) & land == `l', beta
	matrix beta = e(beta)
		replace BSA_b_1622 = _b[year] if land == `l'
		replace BSA_se_1622 = _se[year] if land == `l'
		replace BSA_n_1622 = e(N) if land == `l'
		replace BSA_beta_1622 = beta[1,1] if land == `l'	
		}
	}
}

gen BSec_b_0212 = .
gen BSec_se_0212 = .
gen BSec_n_0212 = .
gen BSec_beta_0212 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSec year if inrange(year,2002,2012) & land == `l', beta
	matrix beta = e(beta)
		replace BSec_b_0212 = _b[year] if land == `l'
		replace BSec_se_0212 = _se[year] if land == `l'
		replace BSec_n_0212 = e(N) if land == `l'
		replace BSec_beta_0212 = beta[1,1] if land == `l'
		}
	}
}

gen BSec_b_1222 = .
gen BSec_se_1222 = .
gen BSec_n_1222 = .
gen BSec_beta_1222 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSec year if inrange(year,2012,2022) & land == `l', beta
	matrix beta = e(beta)
		replace BSec_b_1222  = _b[year] if land == `l'
		replace BSec_se_1222 = _se[year] if land == `l'
		replace BSec_n_1222  = e(N) if land == `l'
		replace BSec_beta_1222 = beta[1,1] if land == `l'
		}
	}
}

gen BSec_b_1622 = .
gen BSec_se_1622 = .
gen BSec_n_1622 = .
gen BSec_beta_1622 = .

foreach l of local countries { 
		quietly {
			capture {
	reg BSec year if inlist(year,2016,2020,2022) & land == `l', beta
	matrix beta = e(beta)
		replace BSec_b_1622 = _b[year] if land == `l'
		replace BSec_se_1622 = _se[year] if land == `l'
		replace BSec_n_1622 = e(N) if land == `l'
		replace BSec_beta_1622 = beta[1,1] if land == `l'
		}
	}
}

save 1_part_one.dta, replace

collapse BS*, by(land)

** Figure 2a
sum BSA_b
tab land if BSA_b < 0
bysort land: tab BSA_n if BSA_b < 0
tab land if BSA_b == .
tab BSA_n
scatter land BSA_b if BSA_n == 11, msym(oh) mlwidth(medium) msize(vhuge) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 10, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 9, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 8, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 7, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 6, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || ///
	scatter land BSA_b if BSA_n == 5, msym(oh) mlwidth(medium) msize(medsmall) mc(edkblue) || /// 
	scatter land BSA_b if BSA_n == 3, msym(oh) mlwidth(medium) msize(small) mc(edkblue) || /// 
	scatter land BSA_b if BSA_n == 2, msym(oh) mlwidth(medium) msize(vsmall) mc(edkblue) || ///
	kdensity BSA_b, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle(" " "Belief System {it:A}" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") xlab(, nolabel noticks) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0024185, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "11" 2 "10" 3 "9" 4 "8" 5 "7" 6 "6" 7 "5" 8 "3" 9 "2" - 10 "Density") stack placement(center) ///
	region(lstyle(none))) name(BSA, replace)
	
sum BSec_b
tab land if BSec_b < 0
bysort land: tab BSec_n if BSec_b < 0
tab land if BSec_b == .
tab BSec_n
scatter land BSec_b if BSec_n == 11, msym(oh) mlwidth(medium) msize(vhuge) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 10, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 9, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 8, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 7, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 6, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || ///
	scatter land BSec_b if BSec_n == 5, msym(oh) mlwidth(medium) msize(medsmall) mc(edkblue) || /// 
	scatter land BSec_b if BSec_n == 3, msym(oh) mlwidth(medium) msize(small) mc(edkblue) || /// 
	scatter land BSec_b if BSec_n == 2, msym(oh) mlwidth(medium) msize(vsmall) mc(edkblue) || ///
	kdensity BSec_b, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle("Economic" "Belief System" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0008412, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "11" 2 "10" 3 "9" 4 "8" 5 "7" 6 "6" 7 "5" 8 "3" 9 "2" - 10 "Density") stack placement(center) ///
	region(lstyle(none))) name(BSec, replace)
	
tab land if BSec_b == . & BSA_b != .
sum BSA_b
local N `r(N)'
grc1leg BSA BSec, xcommon rows(2) title("2002–2022 ({it:n} = `N')")

** Figure 2b
sum BSA_b_0212
tab land if BSA_b_0212 < 0
bysort land: tab BSA_n_0212 if BSA_b_0212 < 0
tab land if BSA_b_0212 == .
tab BSA_n_0212
scatter land BSA_b_0212 if BSA_n_0212 == 6, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSA_b_0212 if BSA_n_0212 == 5, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSA_b_0212 if BSA_n_0212 == 4, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSA_b_0212 if BSA_n_0212 == 3, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSA_b_0212 if BSA_n_0212 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || ///
	kdensity BSA_b_0212, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle(" " "Belief System {it:A}" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") xlab(, nolabel noticks) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0008886, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "6" 2 "5" 3 "4" 4 "3" 5 "2" - 6 "Density") stack placement(center) region(lstyle(none))) name(BSA_0212, replace)
	
sum BSec_b_0212
tab land if BSec_b_0212 < 0
bysort land: tab BSec_n_0212 if BSec_b_0212 < 0
tab land if BSec_b_0212 == .
tab BSec_n_0212
scatter land BSec_b_0212 if BSec_n_0212 == 6, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSec_b_0212 if BSec_n_0212 == 5, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSec_b_0212 if BSec_n_0212 == 4, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSec_b_0212 if BSec_n_0212 == 3, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSec_b_0212 if BSec_n_0212 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || /// 
	kdensity BSec_b_0212, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle("Economic" "Belief System" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0007129, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "6" 2 "5" 3 "4" 4 "3" 5 "2" - 6 "Density") stack placement(center) region(lstyle(none))) name(BSec_0212, replace)
	
tab land if BSec_b_0212 == . & BSA_b_0212 != .
sum BSA_b_0212
local N_0212 `r(N)'
grc1leg BSA_0212 BSec_0212, xcommon rows(2) title("2002–2012 ({it:n} = `N_0212')")

** Figure 2c
sum BSA_b_1222
tab land if BSA_b_1222 < 0
bysort land: tab BSA_n_1222 if BSA_b_1222 < 0
tab land if BSA_b_1222 == .
tab BSA_n_1222
scatter land BSA_b_1222 if BSA_n_1222 == 6, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSA_b_1222 if BSA_n_1222 == 5, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSA_b_1222 if BSA_n_1222 == 4, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSA_b_1222 if BSA_n_1222 == 3, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSA_b_1222 if BSA_n_1222 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || /// 
	kdensity BSA_b_1222, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle(" " "Belief System {it:A}" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") xlab(, nolabel noticks) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0052755, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "6" 2 "5" 3 "4" 4 "3" 5 "2" - 6 "Density") stack placement(center) region(lstyle(none))) name(BSA_1222, replace)
	
sum BSec_b_1222
tab land if BSec_b_1222 < 0
bysort land: tab BSec_n_1222 if BSec_b_1222 < 0
tab land if BSec_b_1222 == .
tab BSec_n_1222
scatter land BSec_b_1222 if BSec_n_1222 == 6, msym(oh) mlwidth(medium) msize(huge) mc(edkblue) || ///
	scatter land BSec_b_1222 if BSec_n_1222 == 5, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSec_b_1222 if BSec_n_1222 == 4, msym(oh) mlwidth(medium) msize(large) mc(edkblue) || ///
	scatter land BSec_b_1222 if BSec_n_1222 == 3, msym(oh) mlwidth(medium) msize(medlarge) mc(edkblue) || ///
	scatter land BSec_b_1222 if BSec_n_1222 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || /// 
	kdensity BSec_b_1222, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle("Economic" "Belief System" " ", size(large)) ysca(off axis(2)) ///
	xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0026771, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "6" 2 "5" 3 "4" 4 "3" 5 "2" - 6 "Density") stack placement(center) region(lstyle(none))) name(BSec_1222, replace)

tab land if BSec_b_1222 == . & BSA_b_1222 != .
sum BSA_b_1222
local N_1222 `r(N)'
grc1leg BSA_1222 BSec_1222, xcommon rows(2) title("2012–2022 ({it:n} = `N_1222')") name(b, replace)

	* t-tests
	sum BSA_n if BSA_b < 0
	sum BSA_n if BSA_b > 0
		gen pos = 1 if BSA_b > 0 & BSA_b != .
		replace pos = 0 if BSA_b < 0
	ttest BSA_n, by(pos)
	tab land BSA_n if pos == 0

	sum BSA_n_0212 if BSA_b_0212 < 0
	sum BSA_n_0212 if BSA_b_0212 > 0
		gen pos_0212 = 1 if BSA_b_0212 > 0 & BSA_b_0212 != .
		replace pos_0212 = 0 if BSA_b_0212 < 0
	ttest BSA_n_0212, by(pos_0212)

	sum BSA_n_1222 if BSA_b_1222 < 0
	sum BSA_n_1222 if BSA_b_1222 > 0
		gen pos_1222 = 1 if BSA_b_1222 > 0 & BSA_b_1222 != .
		replace pos_1222 = 0 if BSA_b_1222 < 0
	ttest BSA_n_1222, by(pos_1222)
	tab land BSA_n_1222 if pos_1222 == 0

** Figure 2d
sum BSA_b_1622
tab land if BSA_b_1622 < 0
tab land if BSec_b_1622 < 0
bysort land: tab BSA_n_1622 if BSA_b_1622 < 0
tab land if BSA_b_1622 == .
tab BSA_n_1622
scatter land BSA_b_1622 if BSA_n_1622 == 3, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || ///
	scatter land BSA_b_1622 if BSA_n_1622 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || ///
	kdensity BSA_b_1622, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle(" " "Belief System {it:A}" " ", size(medlarge)) ysca(off axis(2)) ///
	xsca(off) xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0055877, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "3" 2 "2" - 3 "Density") stack placement(center) region(lstyle(none))) name(BSA_1622, replace)

sum BSec_b_1622
tab land if BSec_b_1622 < 0
bysort land: tab BSec_n_1622 if BSec_b_1622 < 0
tab land if BSec_b_1622 == .
tab BSec_n_1622
scatter land BSec_b_1622 if BSec_n_1622 == 3, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || /// 
	scatter land BSec_b_1622 if BSec_n_1622 == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || /// 
	kdensity BSec_b_1622, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle("Economic" "Belief System" " ", size(medlarge)) ysca(off axis(2)) ///
	xsca(off) xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.0042377, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "3" 2 "2" - 3 "Density") stack placement(center) region(lstyle(none))) name(BSec_1622, replace)
	
sum BSB_b
tab land if BSB_b < 0
bysort land: tab BSB_n if BSB_b < 0
tab land if BSB_b == .
tab BSB_n
scatter land BSB_b if BSB_n == 3, msym(oh) mlwidth(medium) msize(vlarge) mc(edkblue) || /// 
	scatter land BSB_b if BSB_n == 2, msym(oh) mlwidth(medium) msize(medium) mc(edkblue) || /// 
	kdensity BSB_b, fcolor(none) lcolor(edkblue*.5) lwidth(medthin) lp(solid) yaxis(2) ||, ///
	ylab(, nolabel noticks axis(1)) ytitle(" " "Belief System {it:B}" " ", size(medlarge)) ysca(off axis(2)) ///
	xtitle("") ylab(, nogrid) xlab(, labsize(10-pt)) ///
	xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) xline(.008428, lc(cranberry) lwidth(medthin) lp(solid)) ///
	scheme(s1mono) legend(pos(6) rows(1) symplacement(center) size(medsmall) ///
	order(- "{it:t}:" 1 "3" 2 "2" - 3 "Density") stack placement(center) region(lstyle(none))) name(BSB, replace)
	
tab land if BSec_b_1622 == . & BSA_b_1622 != . & BSB_b != .
sum BSA_b_1622
local N_1622 `r(N)'
grc1leg BSA_1622 BSec_1622 BSB, xcommon rows(3) title("2016, 2020, and 2022 ({it:n} = `N_1622')")

tab land if BSA_b_1622 < 0 & BSA_n_1622 == 3
tab land if BSB_b < 0 & BSB_n == 3

// ESTIMATE OVERALL TREND
use 1_part_one.dta, clear

* a) MULTILEVEL MODEL
foreach var of varlist BSA BSB BSec {
	mixed `var' || land:, mle
	estat icc

	mixed `var' year || land:, mle
	eststo `var'
	mixed `var' year || land: year, mle
	eststo `var'_rts
	
	lrtest `var' `var'_rts // LR-Test suggests no better fit with random time slope
}

mixed BSA year || land:, mle var
	eststo ML_BSA
	estat icc
	predict BSA_hat1
	predict BSA_hse1, stdp
	
mixed BSA year if inrange(year,2002,2012) || land:, mle var
	eststo ML_BSA_0212
	estat icc
	predict BSA_hat2
	predict BSA_hse2, stdp

mixed BSA year if inrange(year,2012,2022) || land:, mle var
	eststo ML_BSA_1222
	estat icc
	predict BSA_hat3
	predict BSA_hse3, stdp
	
mixed BSA year if inlist(year,2016,2020,2022) || land:, mle var
	eststo ML_BSA_1622
	estat icc
	predict BSA_hat4
	predict BSA_hse4, stdp
	
mixed BSB year if inlist(year,2016,2020,2022) || land:, mle var
	eststo ML_BSB
	estat icc
	predict BSB_hat
	predict BSB_hse, stdp

mixed BSec year || land:, mle var
	eststo ML_BSec
	estat icc
	predict BSec_hat1
	predict BSec_hse1, stdp

mixed BSec year if inrange(year,2002,2012) || land:, mle var
	eststo ML_BSec_0212
	estat icc
	predict BSec_hat2
	predict BSec_hse2, stdp
	
mixed BSec year if inrange(year,2012,2022) || land:, mle var
	eststo ML_BSec_1222
	estat icc
	predict BSec_hat3
	predict BSec_hse3, stdp
	
mixed BSec year if inlist(year,2016,2020,2022) || land:, mle var
	eststo ML_BSec_1622
	estat icc
	predict BSec_hat4
	predict BSec_hse4, stdp

** Figure 3a
coefplot (ML_BSA, offset(0.1) ms(Oh) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = `" "2002–2022"" " "', angle(vertical) notick labsize(medlarge)) ///
		xsca(extend off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_a, replace)
coefplot (ML_BSA_0212, offset(0.1) ms(Oh) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_0212, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = `" "2002–2012"" " "', angle(vertical) notick labsize(medlarge)) ///
		xsca(extend off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_b, replace)
coefplot (ML_BSA_1222, offset(0.1) ms(Oh) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_1222, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = `" "2012–2022"" " "', angle(vertical) notick labsize(medlarge)) ///
		xsca(extend off) drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_c, replace)
coefplot (ML_BSA_1622, offset(0.25) ms(Oh) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_1622, ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSB, ms(Sh) offset(-0.25) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = `" "2016, 2020,""and 2022" "', angle(vertical) notick labsize(medlarge)) ///
		xsca(extend) drop(_cons) xlab(, nogrid labsize(medlarge)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) plotr(sty(none)) ///
		plotl("Belief System {it:A}" "Economic Belief System" "Belief System {it:B}" ) name(ML_d, replace)

grc1leg ML_a ML_b ML_c ML_d, xcommon cols(1) legendfrom(ML_d) imargin(0 0 0 0) name(ML_AEc, replace)

scatter BSA_hat1

// Robustness checks (for APPENDIX)
coefplot (ML_BSA, offset(0.1) ms(O) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("Multilevel", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_a_app, replace)
coefplot (ML_BSA_0212, offset(0.1) ms(O) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_0212, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("Multilevel", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_b_app, replace)
coefplot (ML_BSA_1222, offset(0.1) ms(O) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_1222, offset(-0.1) ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("Multilevel", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(ML_c_app, replace)
coefplot (ML_BSA_1622, offset(0.2) ms(O) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSec_1622, ms(X) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_BSB, ms(S) offset(-0.2) mc(edkblue) msize(vlarge) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) plotr(sty(none)) ytitle("Multilevel", size(medlarge)) ///
		plotl("Belief System {it:A}" "Economic Belief System" "Belief System {it:B}" ) name(ML_d_app, replace)

bysort land: gen n = _N

foreach var of varlist BSA BSec {
mixed `var' year if inlist(n,1,2) || land:, mle
	eststo ML_`var'2
mixed `var' year if inrange(n,1,3) || land:, mle
	eststo ML_`var'3
mixed `var' year if inrange(n,1,4) || land:, mle
	eststo ML_`var'4
mixed `var' year if inrange(n,1,5) || land:, mle
	eststo ML_`var'5
mixed `var' year if inrange(n,1,6) || land:, mle
	eststo ML_`var'6
mixed `var' year if inrange(n,1,7) || land:, mle
	eststo ML_`var'7
mixed `var' year if inrange(n,1,8) || land:, mle
	eststo ML_`var'8
mixed `var' year if inrange(n,1,9) || land:, mle
	eststo ML_`var'9
mixed `var' year if inrange(n,1,11) || land:, mle
	eststo ML_`var'10
}

foreach var of varlist BSA BSec {
coefplot (ML_`var'2, offset(0.4) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'3, offset(0.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'4, offset(0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'5, offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'6, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'7, offset(-0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'8, offset(-0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'9, offset(-0.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'10, offset(-0.4) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	xsca(extend) drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("Multilevel" " ") name(ML_`var'_rob, replace)
}
		
keep if inrange(year,2002,2012)
bysort land: gen n_0212 = _N

foreach var of varlist BSA BSec {
mixed `var' year if inlist(n_0212,1,2) || land:, mle
	eststo ML_`var'2_0212
mixed `var' year if inrange(n_0212,1,3) || land:, mle
	eststo ML_`var'3_0212
mixed `var' year if inrange(n_0212,1,4) || land:, mle
	eststo ML_`var'4_0212
mixed `var' year if inrange(n_0212,1,5) || land:, mle
	eststo ML_`var'5_0212
mixed `var' year if inrange(n_0212,1,6) || land:, mle
	eststo ML_`var'6_0212
}

foreach var of varlist BSA BSec {
coefplot (ML_`var'2_0212, offset(.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'3_0212, offset(.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'4_0212, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'5_0212, offset(-.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'6_0212, offset(-.3) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("Multilevel" " ") name(ML_`var'_0212_rob, replace)
}

use 1_part_one.dta, clear
keep if inrange(year,2012,2022)
bysort land: gen n_1222 = _N

foreach var of varlist BSA BSec {
mixed `var' year if inlist(n_1222,1,2) || land:, mle
	eststo ML_`var'2_1222
mixed `var' year if inrange(n_1222,1,3) || land:, mle
	eststo ML_`var'3_1222
mixed `var' year if inrange(n_1222,1,4) || land:, mle
	eststo ML_`var'4_1222
mixed `var' year if inrange(n_1222,1,5) || land:, mle
	eststo ML_`var'5_1222
mixed `var' year if inrange(n_1222,1,6) || land:, mle
	eststo ML_`var'6_1222
}

foreach var of varlist BSA BSec {
coefplot (ML_`var'2_1222, offset(.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'3_1222, offset(.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'4_1222, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'5_1222, offset(-.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'6_1222, offset(-.3) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("Multilevel" " ") name(ML_`var'_1222_rob, replace)
}

use 1_part_one.dta, clear
keep if inlist(year,2016,2020,2022)
bysort land: gen n_1622 = _N

foreach var of varlist BSA BSec BSB {
mixed `var' year if inlist(n_1622,1,2) || land:, mle
	eststo ML_`var'2_1622
mixed `var' year if inrange(n_1622,1,3) || land:, mle
	eststo ML_`var'3_1622
}

foreach var of varlist BSA BSec BSB {
coefplot (ML_`var'2_1622, offset(.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(ML_`var'3_1622, offset(-.2) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("Multilevel" " ") name(ML_`var'_1622_rob, replace)
}

* b) GENERALIZED ESTIMATING EQUATIONS
use 1_part_one.dta, clear
collapse BSA BSB BSec, by(land year)

xtset land year
xtgee BSA year, family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA
	estadd beta
	mat GEE_BSA = e(beta)
	matrix list e(beta)
xtgee BSA year if inrange(year,2002,2012), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA_0212
	estadd beta
	mat GEE_BSA_0212 = e(beta)
	matrix list e(beta)
xtgee BSA year if inrange(year,2012,2022), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA_1222
	estadd beta
	mat GEE_BSA_1222 = e(beta)
	matrix list e(beta)
xtgee BSA year if inlist(year,2016,2020,2022), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA_1622
	estadd beta
	mat GEE_BSA_1622 = e(beta)
	matrix list e(beta)
xtgee BSB year if inlist(year,2016,2020,2022), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSB
	estadd beta
	mat GEE_BSB = e(beta)
	matrix list e(beta)
	
xtgee BSec year, family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec
	estadd beta
	mat GEE_BSec = e(beta)
	matrix list e(beta)
xtgee BSec year if inrange(year,2002,2012), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec_0212
	estadd beta
	mat GEE_BSec_0212 = e(beta)
	matrix list e(beta)
xtgee BSec year if inrange(year,2012,2022), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec_1222
	estadd beta
	mat GEE_BSec_1222 = e(beta)
	matrix list e(beta)
xtgee BSec year if inlist(year,2016,2020,2022), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec_1622
	estadd beta
	mat GEE_BSec_1622 = e(beta)
	matrix list e(beta)

coefplot (GEE_BSA, offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_BSec, offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		drop(_cons) xlab(, nogrid labsize(medlarge)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("GEE", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(GEE_a, replace)
coefplot (GEE_BSA_0212, offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_BSec_0212, offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		drop(_cons) xlab(, nogrid labsize(medlarge)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("GEE", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(GEE_b, replace)
coefplot (GEE_BSA_1222, offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_BSec_1222, offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		drop(_cons) xlab(, nogrid labsize(medlarge)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("GEE", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(GEE_c, replace)
coefplot (GEE_BSA_1622, offset(0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_BSB, ms(Sh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_BSec_1622, offset(-0.2) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
		drop(_cons) xlab(, nogrid labsize(medlarge)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(pos(6) rows(1) size(medsmall)) ytitle("GEE", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Belief System {it:B}" "Economic Belief System") name(GEE_d, replace)

// Robustness checks (preparation for online appendix; compiled from line 1283)
bysort land: gen n = _N
xtset land year

**# convergence message -- run models and store estimates line by line
xtgee BSA year if inlist(n,1,2), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA2
xtgee BSA year if inrange(n,1,3), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA3
xtgee BSA year if inrange(n,1,4), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA4
xtgee BSA year if inrange(n,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA5
	
xtgee BSA year if inrange(n,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA6
xtgee BSA year if inrange(n,1,7), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA7
xtgee BSA year if inrange(n,1,8), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA8
xtgee BSA year if inrange(n,1,9), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA9
xtgee BSA year if inrange(n,1,11), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA10

xtgee BSec year if inlist(n,1,2), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec2
xtgee BSec year if inrange(n,1,3), family(gaussian) link(identity) vce(robust)
	**# estimates diverging (correlation > 1) -- run models and store estimates line by line
	eststo GEE_BSec3
xtgee BSec year if inrange(n,1,4), family(gaussian) link(identity) vce(robust)
	**# estimates diverging (correlation > 1) -- run models and store estimates line by line
	eststo GEE_BSec4

xtgee BSec year if inrange(n,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec5
xtgee BSec year if inrange(n,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec6
xtgee BSec year if inrange(n,1,7), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec7
xtgee BSec year if inrange(n,1,8), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec8
xtgee BSec year if inrange(n,1,9), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec9
xtgee BSec year if inrange(n,1,11), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec10

foreach var of varlist BSA BSec {
coefplot (GEE_`var'2, offset(0.4) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'3, offset(0.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'4, offset(0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'5, offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'6, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'7, offset(-0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'8, offset(-0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'9, offset(-0.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'10, offset(-0.4) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	xsca(extend) drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("GEE" " ") name(GEE_`var'_rob, replace)
}
	
keep if inrange(year,2002,2012)
bysort land: gen n_0212 = _N

xtgee BSA year if inlist(n_0212,1,2), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA2_0212
xtgee BSA year if inrange(n_0212,1,3), family(gaussian) link(identity) vce(robust)
**# convergence message -- run further from here
	eststo GEE_BSA3_0212
xtgee BSA year if inrange(n_0212,1,4), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA4_0212
xtgee BSA year if inrange(n_0212,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA5_0212
xtgee BSA year if inrange(n_0212,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA6_0212

xtgee BSec year if inlist(n_0212,1,2), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec2_0212
xtgee BSec year if inrange(n_0212,1,3), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec3_0212
xtgee BSec year if inrange(n_0212,1,4), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec4_0212
xtgee BSec year if inrange(n_0212,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec5_0212
xtgee BSec year if inrange(n_0212,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec6_0212

foreach var of varlist BSA BSec {
coefplot (GEE_`var'2_0212, offset(.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'3_0212, offset(.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'4_0212, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'5_0212, offset(-.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'6_0212, offset(-.3) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("GEE" " ") name(GEE_`var'_0212_rob, replace)
}

use 1_part_one.dta, clear
collapse BSA BSec, by(land year)

xtset land year

keep if inrange(year,2012,2022)
bysort land: gen n_1222 = _N

xtgee BSA year if inlist(n_1222,1,2), family(gaussian) link(identity) vce(robust) 
**# convergence message -- run further from here
	eststo GEE_BSA2_1222
xtgee BSA year if inrange(n_1222,1,3), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA3_1222
xtgee BSA year if inrange(n_1222,1,4), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA4_1222
xtgee BSA year if inrange(n_1222,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA5_1222
xtgee BSA year if inrange(n_1222,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSA6_1222

xtgee BSec year if inlist(n_1222,1,2), family(gaussian) link(identity) vce(robust) 
	eststo GEE_BSec2_1222
xtgee BSec year if inrange(n_1222,1,3), family(gaussian) link(identity) vce(robust)
**# convergence message -- run further from here
	eststo GEE_BSec3_1222
xtgee BSec year if inrange(n_1222,1,4), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec4_1222
xtgee BSec year if inrange(n_1222,1,5), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec5_1222
xtgee BSec year if inrange(n_1222,1,6), family(gaussian) link(identity) vce(robust)
	eststo GEE_BSec6_1222

foreach var of varlist BSA BSec {
coefplot (GEE_`var'2_1222, offset(.3) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'3_1222, offset(.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'4_1222, ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'5_1222, offset(-.15) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'6_1222, offset(-.3) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("GEE" " ") name(GEE_`var'_1222_rob, replace)
}

use 1_part_one.dta, clear
collapse BSA BSB BSec, by(land year)

xtset land year

keep if inlist(year,2016,2020,2022)
bysort land: gen n_1622 = _N

foreach var of varlist BSA BSB BSec {
xtgee `var' year if inlist(n_1622,1,2), family(gaussian) link(identity) vce(robust)
	eststo GEE_`var'2_1622
xtgee `var' year if inrange(n_1622,1,3), family(gaussian) link(identity) vce(robust)
	eststo GEE_`var'3_1622
}

foreach var of varlist BSA BSB BSec {
coefplot (GEE_`var'2_1622, offset(.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
	(GEE_`var'3_1622, offset(-.2) ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(year = " ", notick) ///
	drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		legend(off) plotr(sty(none)) title("GEE" " ") name(GEE_`var'_1622_rob, replace)
}

* c) META-ANALYSIS	
use 1_part_one.dta, clear
collapse *_b* *_se* *_n*, by(land)
	egen allmiss = rowmiss(_all)
	tab land if allmiss == 36 // Albania, Luxembourg, North Macedonia, Romania, Kosovo
	drop if allmiss == 36
	decode land, gen(landstr)
save 2_p1_meta.dta, replace
foreach var of varlist *_se* {
	use 2_p1_meta.dta, replace
	tab land if `var' <= 0
	drop if `var' <= 0
save 2_`var'.dta, replace
}

use 2_BSA_se.dta, clear
meta set BSA_b BSA_se, studylabel(landstr)
meta summarize, random nostudies
matrix BSAmet = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSA_se_0212.dta, clear
meta set BSA_b_0212 BSA_se_0212, studylabel(landstr)
meta summarize, random nostudies
matrix BSAmet0212 = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSA_se_1222.dta, clear
meta set BSA_b_1222 BSA_se_1222, studylabel(landstr)
meta summarize, random nostudies
matrix BSAmet1222 = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSA_se_1622.dta, clear
meta set BSA_b_1622 BSA_se_1622, studylabel(landstr)
meta summarize, random nostudies
matrix BSAmet1622 = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSB_se.dta, clear
meta set BSB_b BSB_se, studylabel(landstr)
meta summarize, random nostudies
matrix BSBmet = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSec_se.dta, clear
meta set BSec_b BSec_se, studylabel(landstr)
meta summarize, random nostudies
matrix BSecmet = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSec_se_0212.dta, clear
meta set BSec_b_0212 BSec_se_0212, studylabel(landstr)
meta summarize, random nostudies
matrix BSecmet0212 = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSec_se_1222.dta, clear
meta set BSec_b_1222 BSec_se_1222, studylabel(landstr)
meta summarize, random nostudies
matrix BSecmet1222 = (r(theta) \ r(ci_lb) \ r(ci_ub))

use 2_BSec_se_1622.dta, clear
meta set BSec_b_1622 BSec_se_1622, studylabel(landstr)
meta summarize, random nostudies
matrix BSecmet1622 = (r(theta) \ r(ci_lb) \ r(ci_ub))

coefplot (mat(BSAmet[1,.]), ci((BSAmet[2,.] BSAmet[3,.])) ///
			offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmet[1,.]), ci((BSecmet[2,.] BSecmet[3,.])) ///
			offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(pos(6) rows(1)) ytitle("Meta-analysis", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(MET_a, replace)
		
coefplot (mat(BSAmet0212[1,.]), ci((BSAmet0212[2,.] BSAmet0212[3,.])) ///
			offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmet0212[1,.]), ci((BSecmet0212[2,.] BSecmet0212[3,.])) ///
			offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(pos(6) rows(1)) ytitle("Meta-analysis", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(MET_b, replace)
		
coefplot (mat(BSAmet1222[1,.]), ci((BSAmet1222[2,.] BSAmet1222[3,.])) ///
			offset(0.1) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmet1222[1,.]), ci((BSecmet1222[2,.] BSecmet1222[3,.])) ///
			offset(-0.1) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(pos(6) rows(1)) ytitle("Meta-analysis", size(medlarge)) ///
		plotr(sty(none)) plotl("Belief System {it:A}" "Economic Belief System") name(MET_c, replace)

coefplot (mat(BSAmet1622[1,.]), ci((BSAmet1622[2,.] BSAmet1622[3,.])) ///
			offset(0.2) ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSBmet[1,.]), ci((BSBmet[2,.] BSBmet[3,.])) ///
			ms(Sh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmet1622[1,.]), ci((BSecmet1622[2,.] BSecmet1622[3,.])) ///
			offset(-0.2) ms(X) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", notick) ///
		xsca(off) drop(_cons) xlab(, nogrid) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(pos(6) rows(1)) ytitle("Meta-analysis", size(medlarge)) ///
		plotl("Belief System {it:A}" "Belief System {it:B}" "Economic Belief System") name(MET_d, replace)

// Robustness checks (preparation for online appendix; compiled from line 1283)
use 2_BSA_se.dta, clear
meta set BSA_b BSA_se, studylabel(landstr)

tab BSA_n
gen sg1 = 1 if BSA_n == 3
gen sg2 = 1 if inlist(BSA_n,3,5)
gen sg3 = 1 if inlist(BSA_n,3,6)
gen sg4 = 1 if inrange(BSA_n,3,7)
gen sg5 = 1 if inrange(BSA_n,3,8)
gen sg6 = 1 if inrange(BSA_n,3,9)
gen sg7 = 1 if inrange(BSA_n,3,10)
gen sg8 = 1 if inrange(BSA_n,3,11)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSAmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSAmetsg1[1,.]), ci((BSAmetsg1[2,.] BSAmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg2[1,.]), ci((BSAmetsg2[2,.] BSAmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg3[1,.]), ci((BSAmetsg3[2,.] BSAmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg4[1,.]), ci((BSAmetsg4[2,.] BSAmetsg4[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg5[1,.]), ci((BSAmetsg5[2,.] BSAmetsg5[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg6[1,.]), ci((BSAmetsg6[2,.] BSAmetsg6[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg7[1,.]), ci((BSAmetsg7[2,.] BSAmetsg7[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg8[1,.]), ci((BSAmetsg8[2,.] BSAmetsg8[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSA_rob, replace)

use 2_BSA_se_0212.dta, clear
meta set BSA_b_0212 BSA_se_0212, studylabel(landstr)

tab BSA_n_0212
gen sg1 = 1 if BSA_n_0212 == 3
gen sg2 = 1 if inlist(BSA_n_0212,3,4)
gen sg3 = 1 if inlist(BSA_n_0212,3,5)
gen sg4 = 1 if inrange(BSA_n_0212,3,6)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSAmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSAmetsg1[1,.]), ci((BSAmetsg1[2,.] BSAmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg2[1,.]), ci((BSAmetsg2[2,.] BSAmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg3[1,.]), ci((BSAmetsg3[2,.] BSAmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg4[1,.]), ci((BSAmetsg4[2,.] BSAmetsg4[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSA_0212_rob, replace)

use 2_BSA_se_1222.dta, clear
meta set BSA_b_1222 BSA_se_1222, studylabel(landstr)

tab BSA_n_1222
gen sg1 = 1 if BSA_n_1222 == 3
gen sg2 = 1 if inlist(BSA_n_1222,3,4)
gen sg3 = 1 if inrange(BSA_n_1222,3,5)
gen sg4 = 1 if inrange(BSA_n_1222,3,6)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSAmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSAmetsg1[1,.]), ci((BSAmetsg1[2,.] BSAmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg2[1,.]), ci((BSAmetsg2[2,.] BSAmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg3[1,.]), ci((BSAmetsg3[2,.] BSAmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSAmetsg4[1,.]), ci((BSAmetsg4[2,.] BSAmetsg4[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSA_1222_rob, replace)

use 2_BSec_se.dta, clear
meta set BSec_b BSec_se, studylabel(landstr)

tab BSec_n
gen sg1 = 1 if BSec_n == 3
gen sg2 = 1 if inlist(BSec_n,3,5)
gen sg3 = 1 if inlist(BSec_n,3,6)
gen sg4 = 1 if inrange(BSec_n,3,7)
gen sg5 = 1 if inrange(BSec_n,3,8)
gen sg6 = 1 if inrange(BSec_n,3,9)
gen sg7 = 1 if inrange(BSec_n,3,10)
gen sg8 = 1 if inrange(BSec_n,3,11)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSecmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSecmetsg1[1,.]), ci((BSecmetsg1[2,.] BSecmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg2[1,.]), ci((BSecmetsg2[2,.] BSecmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg3[1,.]), ci((BSecmetsg3[2,.] BSecmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg4[1,.]), ci((BSecmetsg4[2,.] BSecmetsg4[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg5[1,.]), ci((BSecmetsg5[2,.] BSecmetsg5[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg6[1,.]), ci((BSecmetsg6[2,.] BSecmetsg6[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg7[1,.]), ci((BSecmetsg7[2,.] BSecmetsg7[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg8[1,.]), ci((BSecmetsg8[2,.] BSecmetsg8[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSec_rob, replace)

use 2_BSec_se_0212.dta, clear
meta set BSec_b_0212 BSec_se_0212, studylabel(landstr)

tab BSec_n_0212
gen sg1 = 1 if BSec_n_0212 == 3
gen sg2 = 1 if inlist(BSec_n_0212,3,4)
gen sg3 = 1 if inlist(BSec_n_0212,3,5)
gen sg4 = 1 if inrange(BSec_n_0212,3,6)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSecmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSecmetsg1[1,.]), ci((BSecmetsg1[2,.] BSecmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg2[1,.]), ci((BSecmetsg2[2,.] BSecmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg3[1,.]), ci((BSecmetsg3[2,.] BSecmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg4[1,.]), ci((BSecmetsg4[2,.] BSecmetsg4[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSec_0212_rob, replace)

use 2_BSec_se_1222.dta, clear
meta set BSec_b_1222 BSec_se_1222, studylabel(landstr)

tab BSec_n_1222
gen sg1 = 1 if BSec_n_1222 == 3
gen sg2 = 1 if inlist(BSec_n_1222,3,4)
gen sg3 = 1 if inrange(BSec_n_1222,3,5)
gen sg4 = 1 if inrange(BSec_n_1222,3,6)

foreach var of varlist sg* {
	meta summarize if `var' == 1, random nostudies
	matrix BSecmet`var' = (r(theta) \ r(ci_lb) \ r(ci_ub))
}

coefplot (mat(BSecmetsg1[1,.]), ci((BSecmetsg1[2,.] BSecmetsg1[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg2[1,.]), ci((BSecmetsg2[2,.] BSecmetsg2[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg3[1,.]), ci((BSecmetsg3[2,.] BSecmetsg3[3,.])) ///
			ms(Oh) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))) ///
		(mat(BSecmetsg4[1,.]), ci((BSecmetsg4[2,.] BSecmetsg4[3,.])) ///
			ms(O) mc(edkblue) msize(large) ciopts(recast(rcap) lc(edkblue*.5))), ///
	coeflabels(c1 = " ", angle(vertical) notick labsize(medlarge)) ///
		drop(_cons) xlab(, nogrid labsize(medium)) xline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
		subtitle(, bcolor(white) size(medlarge) justification(left)) legend(off) ///
		plotr(sty(none)) title("Meta-analysis" " ") name(MET_BSec_1222_rob, replace)

	// APPENDIX
	// Robustness checks
	** Figure A2a
	grc1leg ML_a_app MET_a GEE_a, rows(3) xcommon imargin(0 0 0 0) legendfrom(GEE_a) title("2002–2022" " ")

	** Figure A2b
	gr combine ML_BSA_rob MET_BSA_rob GEE_BSA_rob, rows(1) xcommon

	** Figure A2c
	gr combine ML_BSec_rob MET_BSec_rob GEE_BSec_rob, rows(1) xcommon

	** Figure A3a
	grc1leg ML_b_app MET_b GEE_b, rows(3) xcommon imargin(0 0 0 0) legendfrom(GEE_b) title("2002–2012" " ")

	** Figure A3b
	gr combine ML_BSA_0212_rob MET_BSA_0212_rob GEE_BSA_0212_rob, rows(1) xcommon

	** Figure A3c
	gr combine ML_BSec_0212_rob MET_BSec_0212_rob GEE_BSec_0212_rob, rows(1) xcommon

	** Figure A4a
	grc1leg ML_c_app MET_c GEE_c, rows(3) xcommon imargin(0 0 0 0) legendfrom(GEE_c) title("2012–2022" " ")

	** Figure A4b
	gr combine ML_BSA_1222_rob MET_BSA_1222_rob GEE_BSA_1222_rob, rows(1) xcommon

	** Figure A4c
	gr combine ML_BSec_1222_rob MET_BSec_1222_rob GEE_BSec_1222_rob, rows(1) xcommon

	** Figure A5a
	grc1leg ML_d_app MET_d GEE_d, rows(3) xcommon imargin(0 0 0 0) legendfrom(GEE_d) title("2016, 2020, and 2022" " ")

	** Figure A5b
	gr combine ML_BSA_1622_rob GEE_BSA_1622_rob, rows(1) xcommon

	** Figure A5c
	gr combine ML_BSec_1622_rob GEE_BSec_1622_rob, rows(1) xcommon

	** Figure A5d
	gr combine ML_BSB_1622_rob GEE_BSB_1622_rob, rows(1) xcommon